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PHASE  SPACE  METHODS  AND  PATH  INTEGRATION: 

A  MICROSCOPIC  APPROACH  TO  OIRECT  AND  INVERSE  WAVE  PROPAGATION 


Louis  Fishman 
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The  Catholic  University  of  America 
Washington,  D.C.  20064  USA 


ABSTRACT.  This  project  focuses  on  the  development  of  new,  multldlmen- 
slonal  algorithms  for  direct  acoustic  propagation  and  generalized  acoustic 
tomography  at  the  level  of  the  scalar  Helmholtz  equation.  The  general  aim 
Is  the  continued  detailed  development  of  the  Ideas  originally  outlined  sev¬ 
eral  years  ago.  Phase  space,  or  "microscopic,"  methods  and  path  (function¬ 
al)  Integral  representations  provide  the  appropriate  framework  to  extend 
homogeneous  Fourier  methods  to  Inhomogeneous  environments.  The  path 
Integrals  furnish  the  principal  representation  of  the  Helmholtz  propagator 
and,  subsequently,  through  direct  computation,  the  basis  for  the  direct 
numerical  algorithms.  There  are  two  complementary  approaches  to  the 
analysis  and  computation  of  the  n-dlmenslonal  Helmholtz  propagator.  The 
first  Is  essentially  a  factorlzatlon/parabollc-based  (one-way)  phase  space 
path  Integration/ Invariant  Imbedding  approach.  This  results  In  a  marching 
algorithm  which  generalizes  the  Tappert/Hardln  split-step  FFT  algorithm  for 
one-way  wave  propagation,  a  nonperturbatlve  Incorporation  of  backscatter 
effects  which  generalizes  Kennett's  algorithm  In  reflection  seismology  for 
two-way  wave  propagation,  and  the  basis  for  the  formulation  and  solution  of 
corresponding  arbitrary-dimensional  nonlinear  Inverse  problems.  The 
numerical  algorithms  based  on  these  modem,  "microscopic"  methods  directly 
compute  pseudo-differential  and  Fourier  Integral  operators.  Incorporate 
phase  space  filtering,  and  are  Ideally  suited  for  computers  which  provide 
either  a  vector  or  a  parallel  pipe  type  of  operation.  Extensive  testing 
has,  so  far,  been  very  promising.  While  the  first  approach  starts  from  a 
transversely  Inhomogeneous  formulation  and,  subsequently,  builds  In 
backscatter  effects,  the  second  approach  constructs  elliptic-based 
(two-t(ay)  P*th  Integral  representations  of  the  propagator  for  general 
range- dependent  environments  from  the  outset.  A  particular  approximate 
path  Integral  construction  (Feynman/Garrod)  results  In  a  true  path 
functional,  suggesting  the  underlying  stochastic  foundations  of  the 
Helmholtz  equation.  It  appears  to  be  a  viable  computational  approximation 
for  a  useful  range  of  propagation  experiments  and  can  be  numerically 
evaluated  by  standard  Monte  Carlo  (statistical)  methods.  A  more  detailed 
examination  and  approximate  construction  of  the  underlying  stochastic 
process  would  provide  for  both  more  accurate  and  widely  applicable  path 
Integral  representations  and  direct  numerical  simulation  techniques. 

I.  INTRODUCTION.  Direct  wave  propagation  modeling  plays  a 
significant  role  in  such  fields  as  underwater  communication,  radio 
transmission  through  the  atmosphere,  laser  propagation,  and  earthquake 
prediction.  Likewise,  the  corresponding  inverse  problems  are  at  the  heart 
of  such  areas  as  submarine  detection,  CAT  scan  technology,  soft- tissue 
diffraction  tomography,  the  mapping  of  the  Interior  earth,  and  oil 


exploration.  In  all  of  these  and  maqy  other  examples,  relatively  fast  and 
accurate  numerical  algorithms  are  necessary. 


The  analysis  and  fast,  accurate  numerical  computation  of  the  wave 
equations  of  classical  physics  are  often  quite  difficult  for  rapidly 
changing,  multidimensional  environments  extending  over  maqy  wavelengths. 

For  the  most  part,  classical,  "macroscopic"  methods  have  resulted  In  direct 
wave  field  approximations  (perturbation  theory,  ray-theory  asymptotics, 
modal  analysis,  hybrid  ray-mode  methods),  derivations  of  approximate  wave 
equations  (scaling  analysis,  field  splitting  techniques,  formal  operator 
expansions),  and  discrete  numerical  approximations  (finite  differences, 
finite  elements,  spectral  methods).  In  the  last  several  decades,  however, 
mathematicians  studying  linear  partial  differential  equations  have 
developed.  In  the  language  of  physicists,  a  sophisticated,  "microscopic" 
phase  space  analysis.  In  conjunction  with  the  global  functional  Integral 
techniques  pioneered  by  Wiener  (Brownian  motion)  and  Feynman  (quantum 
mechanics),  and  so  successfully  applied  todey  In  quantum  field  theory  and 
statistical  physics,  the  n-dlmenslonal  classical  physics  propagators  can  be 
both  represented  explicitly  and  computed  directly.  The  phase  space,  or 
"microscopic,"  methods  and  path  (functional)  Integral  representations 
provide  the  appropriate  framework  to  extend  homogeneous  Fourier  methods  to 
Inhomogeneous  environments.  In  addition  to  suggesting  the  basis  for  the 
formulation  and  solution  of  corresponding  arbitrary-dimensional  nonlinear 
Inverse  problems.  Moreover,  It  Is  In  phase  space,  rather  than  In 
configuration  space,  that,  from  a  mathematical  perspective,  the  Interesting 
geometry  takes  place. 

II.  PHASE  SPACE  ANO  PATH  INTEGRAL  CONSTRUCTIONS.  For  the 
n-dlmenslonal  scalar  Helmholtz  equation,  tnere  are  two  complementary  ap¬ 
proaches  to  this  snslysls  and  computation,  as  Illustrated  In  Figure  1.  The 
first  Is  essentially  a  factorization/path  Integration/Invariant  Imbedding 
approach.  For  transversely  Inhomogeneous  environments.  Implying  medium 
homogeneity  with  respect  to  a  single  distinguished  direction,  the  n- 
dlmenslonal  Helmholtz  equation  can  be  exactly  factored  Into  separate, 
physical  forward  and  backward,  one-way  wave  equations,  following  from 
spectral  analysis  [1-5].  The  forward  evolution  (one-wqy)  equation 

(1/K)dx  *+(*.*t)  ♦  («2(xt)  ♦  (1/R2tft2)l/2  #*(x,xt)  -  0  ,  (1) 

where  K(x)  Is  the  refractive  Index  field  and  E  Is  a  reference  wave  number. 

Is  the  formally  exact  wave  equation  for  propagation  In  a  transversely  In¬ 
homogeneous  half-space  supplemented  with  appropriate  outgoing  wave  radiation 
and  Initial-value  conditions.  While  functions  of  a  finite  set  of  commuting 
self-adjoint  operators  can  be  defined  through  spectral  theory,  functions  of 
noncommuting  operators  arc  represented  by  pseudo-differential  operators 
[2,5].  The  formal  wave  equation  (1)  Is  now  written  explicitly  as  a  Weyl 
pseudo-differential  equation  In  the  form 


(1/E)3x  **(x,xt)  ♦  (E/2w)n_1 


d*idP t 


°l{£t,(-t  *  *t)/2)  exp(il£t-(*t  -  *{))♦*(*, 


**)  -  0. 


(2) 


WAY  ALGORITHMS 


In  Eq.(2),  the  symbol  %(£,q)  associated  with  the  square  root  Helmholtz  oper¬ 
ator  B  »  (k2(<j)  +  (1/t  satisfies  the  Weyl  composition  equation 


K2{£)  -  £2 


(lc/w)2n“2 


dtdxc^dz  Qg(  t+£,  x+£) 


*°B^+£*  l+£)  exp(21t(x‘^  -  t*z))  (3) 

2 

with  Qg2(j),£)  the  symbol  associated  with  the  square  of  B,  B  = 

(K2(q)  ♦  (l/lc2>V  2)  [2,3,5].  The  generalized  Fourier  construction  procedure 

for  the  square  root  Helmholtz  operator  can  be  summarized  plctorlally  by  the 
following  correspondence  diagram 

B2  <=»  Qj.2 

t  r 

B  «=>  PB 

where  the  arrows  symbolize  the  one-  and  two-way  mappings  between  the  appro¬ 
priate  quantities. 

Exact  solutions  of  the  Weyl  composition  equation  (3)  can  be  constructed 
In  several  cases  [6].  For  example,  the  symbol  Qg(p,q)  for  the  two-dimen¬ 
sional  (n  «  2)  quadratic  medium,  K2{q)  ■  Kq  +  w2q2.  Is  given  by  [6] 


Qg(p,q)  *  -(exp(1w/4)t1/2/w1/2)  J  dt  exp(1(Yt  ♦  Xtanht)) 

•t"1^2  (lYsecht  ♦  1Xsech3t  -  (sechtMtanht))  (4) 

with  X  •  (l/t)(w2q2  -  p2),  Y  »  K^/t,  and  t  ■  w/(c.  Consistent  with  taking 

the  square  root  of  the  Indefinite  Helmholtz  operator,  the  corresponding 
symbols,  generally,  have  both  real  and  Imaginary  parts  characterized  by 
oscillatory  behavior  [4,6],  as  Illustrated  In  Figure  2.  Nonuniform  and 
uniform  perturbation  solutions  corresponding  to  definite  physical  limits 
(frequency,  propagation  angle,  field  strength,  field  gradient)  recover 
several  known  approximate  wave  theories  (ordinary  parabolic,  range- 
refraction  parabolic,  Grandvull lemln-extended  parabolic,  half-space  Born, 
Thomson-Chapman ,  rational  linear)  and  systematically  lead  to  several  new 
full-wave,  wide-angle  approximations  [2-4,6]. 

The  exact  pseudo- differential  evolution  equation  (2)  and.  In  general, 
the  wide-angle  extended  parabolic  approximate  equations  derived  from  the 
ana'ysls  of  the  composition  equation  [2-4,6]  are  singular  Integro- 
dlfferentlal  wave  equations.  Solution  representations  for  such  pseudo- 
differential  equations  can  be  directly  expressed  In  terms  of  Infinite- 
dimensional  functional,  or  path,  Integrals  [7,8],  following  from  the  Harkov 
property  of  the  propagator.  In  an  operator  notation,  then, 
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exp(lKBx)  *  11m  TT  exp(ikBAx-)  (5) 

N — >-  j*l  J 

where  AXj  a  x/N,  symbolically  representing  the  propagator  in  terms  of  the 

Infinitesimal  propagator.  As  the  operator  symbol  is  not  simply  quadratic  in 
g,  the  configuration  space  Feynman  path  integral  formulation  is  not  appro¬ 
priate,  necessitating  the  more  general  phase  space  construction  [4,7].  This 
results  in  a  parabolic-based  (one-way)  Hamiltonian  phase  space  path  integral 
representation  of  the  propagator  in  the  form  [3,7] 

,  r  n-i  n  _ , 

G  (x.xJO.x!)  *  lim  /  TT  dx-.  TT  ( C/2ir)  do. 

-t  N - >-  /  j=l  JZ  j=l 

J  R(n-1) (2N-1 ) 

•expfit  J-  (£jt-(xjt  -  Xj.lt)  »  (x/N)  H(£Jt,xjt,Xj.lt)))  (6) 


where 


Fig.  2.  The  real  ( - )  and  imaginary  ( - )  parts  of  the  n  =  2 

quadratic  medium  symbol  as  a  function  of  X  for  '  *  l. 


H (£,<!'  '  )  «  (t/2w) 


dsdt  F(£' -£",£) 

2n-2 

’+£'  )/Z)  ~  t)  exp(iks*t).  (7) 
In  Eq.(7),  F(u,v)  and  hg(£,£)  are  related  to  the  operator  symbol  Og(£,£)  by 

Qg(utv)  *  F(u,v)  hB(u,v)  (8) 

where  fig  and  hg  are  the  corresponding  Fourier  transforms  [2,3,7]. 

The  nonuniqueness  of  the  lattice-approximation  path  integral  represen¬ 
tation  is  readily  understood  In  terms  of  different  discretizations,  or  quad¬ 
ratures,  of  the  symbolic  functional  Integral  and  corresponds  to  the  repre¬ 
sentation  of  a  given  (fixed)  operator  by  different  operator-ordering,  or 
pseudo-differential  operator,  schemes  [2, 3, 7, 8].  More  fundamentally,  in 
analogy  with  the  Schrfidinger  equation  for  particle  motion  on  a  Rlemannian 
space  and  the  thermodynamic  (Fokker-Planck)  equation  for  particle  diffusion, 
the  algorithmic  Helmholtz  path  integral  construction  reflects  the  stochastic 
nature  of  the  Integration  [4,9].  Further,  both  the  macroscopic  and  micro¬ 
scopic  (Infinitesimal)  half-space  propagators  can  be  formally  expressed  as 
Fourier  Integral  operators  with  complex  phase  [4].  The  phase  space  path 
Integral,  thus,  represents  the  macroscopic  Fourier  Integral  operator  in 
terms  of  the  N-fold  application  of  the  microscopic,  or  infinitesimal, 

Fourier  Integral  operator  In  a  manner  which  can  be  related  to  the  global 
geometrical -optics  construction  of  the  macroscopic  operator  [4,5]. 

The  path  integral  formulation  Interprets  the  wave  theory  in  terms  of  an 
Infinitesimal  propagator  summed  over  all  phase  space  paths.  For  the  Helm¬ 
holtz  theory,  the  exact  Infinitesimal  propagator  is  not,  in  general,  given 
by  the  locally  homogeneous  medium  propagator,  as  In  the  ordinary  parabolic 
(Schrbdlnger)  propagator  construction  [8].  The  approximate  extended  para¬ 
bolic  wave  theories  then  correspond  to  approximate  Infinitesimal  propagators 
summed  over  the  complete  phase  space.  In  retaining  the  "sum  over  all 
paths,"  diffraction,  or  full -wave,  effects  are  i ncorporated . 

For  weakly  range-dependent  environments,  range  variability  can  be,  at 
first,  accommodated  at  the  level  of  range  updating,  as  in  the  case  of  the 
parabolic  path  Integral  [1,8],  For  reflection/transmission  from  a  planar 
Interface  separating  two  (different)  transversely  Inhomogeneous  acoustic 
half-spaces,  the  concept  of  reflection  and  transmission  amplitudes  general¬ 
izes  to  reflection  (r)  and  transmission  (t)  operators.  The  reflection  and 
transmission  operators,  which,  when  applied  to  the  incident  wave  field  at 
the  interface,  produce  the  Initial  values  of  the  reflected  and  transmitted 
wave  fields,  are  defined  within  the  Weyl  pseudo-differential  operator 
framework  and  are  explicitly  determined  by  enforcing  the  well-known 
Interface  continuity  conditions.  The  main  result  [10]  Is  a  composition 
equation  of  the  form 


dtdxdydz  x+£)  ^ 

-4 
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GBLt£*S)  _  °br(£»S)  *  (R/w)2n-2 


^  *.«  V  */  >  V  ■✓  *.*'*>'* 

. 


O^tt^,  ®rCjr*£.  *♦£)  exp(2lE(x«£  -  t*z))  (9) 

for  the  reflection  operator  symbol  flr(£,^)  and  an  analogous  equation  for  the 

transmission  operator  symbol  &{(£»£)•  The  Inclusion  of  a  planar  transition 

region  of  arbitrary  length  and  Inhomogeneity  can  be  accomplished  by  factor¬ 
ization  Methods  In  conjunction  with  Invariant  iMbeddlng  [4,11].  Invariant 
iMbeddlng  constructs  the  initial-value  system  for  the  reflection  and  trans¬ 
mission  operators  associated  with  the  transition  region,  transforming  the 
Helmholtz  boundary -value  problem  Into  an  initial -value  problem.  A  dis¬ 
cretized  formulation  [11]  provides  the  extension  of  Kennett's  method  [4,11] 
in  reflection  seismology.  The  resultant  forward  and  backward  wave  fields 
propagating  In  the  transversely  Inhomogeneous  half- spaces  are  represented  by 
the  one-way  path  Integrals,  while,  within  the  transition  region,  a  formal 
path  integral  representation  of  the  propagator  can  be  expressed  as  a  product 
integral  [8].  This  takes  the  form  [4] 

u 

P 

G  *  /  exp( IkH(s)ds)  *  11m  TT  exp(  i(cH(s  - ) As  - )  (10) 

J  “  N - >-  j=l  =  J  J 

where  Sj  *  a  +  (j-l/2)ASj,  ASj  *  (x-a)/N,  a  denotes  the  transition  region 

boundary,  H  is  the  appropriate  first-order  Helmholtz  equation  matrix 

operator  [2,4],  and  with  the  product  of  exponential  factors  ordered  from 
right  (lower  j)  to  left  (higher  j)  reflecting  the  noncommutativity  of  the 
matrix  operator  H  at  different  x.  While  product  Integration-based  path 

integral  constructions  have  been  applied  to  the  problems  of  nonrelativistic 
electron  spin  and  the  Dirac  equation,  such  Infinite  products  of  matrices 
are,  generally,  only  tractable  In  simple  limiting  cases  [4,8]. 


Rather  than  starting  from  a  transversely  Inhomogeneous  formulation  and, 
subsequently,  building  In  backscatter  effects,  the  generalization  of  Fourier 
methods  to  arbitrary  Inhomogeneous  environments  and  the  construction  of  a 
dynamical  basis  for  the  Helmholtz  equation  can  proceed.  In  the  second  ap¬ 
proach,  from  the  construction  of  truly  global  configuration  space  path  inte¬ 
grals,  which  attempt  to  generalize,  for  example,  the  homogeneous  half-space 
result  [3,7] 


G  (x.x^lO.x^)  *  lim 


N-l 

TT  (ixxN 


(n-l )N/2 


j=l 


-jt 


R(n-1 )(N-1) 


(tV2*V 


1)N+1 


( (n-l )N+1 )/2  u(l) 


H( (n-l )N+1 )/2(tK(A  n-l )N+1 } } 


(11) 


where 
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(12) 


^  (n-1  )N*1  ‘  (-Jt  *  -J-lt)2  *  *Z)1/2 

and  H^(|)  Is  the  Hankel  function.  These  elliptic-based  (two-way)  con¬ 
structions,  originating  from  the  Fourier  transform  relationship  between  the 
Helmholtz  and  Schrfldinger  (parabolic)  propagators,  result  In  the  approximate 
Feynman/Garrod  path  Integral  C3, 7 

o  P  N-l  N  exp(  llcSy ) 

G(xlx’)  s  (-1/2^)  11m  /  IT  dx.  TT  (C/2«)ndp.  - —  (13) 

H - >-  J  j=l  j  =  l  J  (1/2  -  r) 

J  Rn(2N-l) 

where 

SN  =  »4> 

corresponds  to  an  appropriate  discretized  action  and 


E  =  (1/N )  (pj/2  +  VUj))  (15) 


plays  a  role  analogous  to  an  average  energy  with  the  Identification  V(x)  2 
2  ~ 
(-1/2HK  (x)  -  1).  For  a  transversely  Inhomogeneous  half-space,  partial 

Integration  of  Eq . ( 13 )  In  conjunction  with  the  reflection  principle  (or 
method  of  images)  results  In  [3,7] 


G+(x,xJ  0,x‘ )  11m 

N - >- 


( n-1 ) (2N-1 ) 


N-l  N 

TT  dxit  TT  (K/2w)n"1d£,t 

j*l  j=l 


•exp(lE(SN  +  21/2x( 1/2  -  E)1/2))  (16) 

with  SN  and  E  taking  on  their  appropriate  forms  In  one- lower  dimension. 

Formally  reducing  both  the  full-  and  transversely  Inhomogeneous 
half-space  phase  space  Feynman/Garrod  path  integrals  to  configuration  space 
path  Integrals  [7]  establishes  the  path  functional  character  of  the 
representation.  Moreover,  the  approximate  Feynman/Garrod  path  Integral  Is 
exact  In  the  homogeneous  medium  limit.  Incorporates  significant  backscatter 
Information,  and  contains  both  the  geometrical  (ray)  acoustic  and  ordinary 
parabolic  approximations.  This  configuration  space  formulation  for  the 
two-way  problem.  Initially  based  on  a  variational  principle  and  phase  space 
constructions,  seeks  to  express  the  propagator  in  terms  of  a  phase 
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functional  evaluated  over  an  appropriate  path  space,  as  symbolically 
expressed  In  the  Feynman/DeWitt-Morette  representation  [3,7,9].  This  takes 
the  form 


G(xlx')  =  (-l/2lc2)  /  D(?)  exp(  itfl(?  ) ) 


where 


lid? II  (1  -  2V(?))J 


is  the  analog  of  the  action  associated  with  a  "free  particle"  on  a  space 
w i th  the  metric 

dl2  =  (1  -  2V (?) ) II  d?l|  2  ( 

and  where  E  represents  the  space  of  paths  from  x'  to  x  such  that 


1/2  =  (1/r)  /  dt  ((1/2)11  d?(t)/dt!r  +  V(?(t))) 


with  the  constraints 
1(0)  =  x-  . 

?(r)  =  x  .  (21) 

The  dynamical  basis  of  the  Helmholtz  equation  can,  thus,  be  viewed  in  terms 
of  a  stochastic  process  embodying  fixed  "average  energy"  paths,  or, 
alternatively,  in  terms  of  "free  particle"  motion  [3,7,9]. 

III.  COMPUTATIONAL  ALGORITHMS.  Direct  integration  of  the  one-way 
phase  space  path  integral  provides  the  computational  basis  for  the  pseudo- 
differential  wave  equation  (2).  Choosing  the  standard  ordering,  F(u,v)  * 
exp(-ilcu*v/2) ,  in  Eqs.  (6),  (7),  and  (8)  results  in  a  numerically  more 
efficient  post-point  marching  algorithm  in  the  form 

4*(x*Ax,xt)^  f  dpt  exp{Uc£toxt)  (exp(1kAxhB(£t,xt))  <t>  +(x,£t))  (22) 


where  is  the  Fourier-transformed  wave  field  and 


hB(£t»*t)  =  (*/»)n_1  /  dsdt  Ogfs.t)  exp(-21R(*t  "  -*‘*£t  ‘  *23* 


This  marching  algorithm  provides  the  generalization  of  the  Tappert/Hardln 
split-step  FFT  algorithm  [1]  to  the  full  one-wty  (factored  Helmholtz)  wave 
equation.  For  a  two-dimensional  model  ocean/bottom  propagation  environment 
with  a  perfectly  reflecting  ocean  surface,  the  Fourier  transform  of  the  wave 
field  In  Eq . ( 22 )  is  replaced  by  a  discrete  fast  sine  transform  and  the  in¬ 
verse  transform  is  evaluated  by  a  rectangular  rule  integration,  enabling  the 
propagated  wave  field  to  be  expressed  in  the  matrix  form 

^+(x+4x,zn)  =  XZ  (24) 

m 

,  + 

for  each  depth  point  zn.  In  Eq.(24),  <p  and  9  are  column  vectors  and  the 
matrix  A  is  defined  by  its  matrix  elements 

Anm  *1sin<EVn  +  EAxhB(‘V2n)  1  1251 

where  hg  and  hg  are  the  even  and  odd  parts  with  respect  to  p  of  hg(p,z)  in 
Eq . ( 23 )  and  Tj  is  an  appropriate  transform  normalization  constant  [1,4,12]. 

The  principal  idea  underlying  the  practical  implementation  of  the  phase 
space  marching  algorithm  is  the  construction  of  a  small  number  of  approxi¬ 
mate  operator  symbols,  which,  when  taken  together,  allow  for  wave  field 
computations  over  a  very  wide  range  of  model  environments  and  propagation 
parameters.  In  conjunction  with  a  study  of  exactly  soluble  cases  of  the 
Weyl  composition  equation  [6],  high-frequency,  real  Weyl  high-frequency, 
uniform  high-frequency,  and  1 ow- f requency  approximate  symbols  have  been  con¬ 
structed  [2-4,6].  Of  particular  significance  is  the  fact  that  the  manner  of 
marching  the  radiation  field  is  independent  of  the  medium  and  any  approxi¬ 
mation  to  the  square  root  Helmholtz  operator,  resulting  in  a  modular  code 
architecture  and  highly  versatile  propagation  program.  Moreover,  the  propa¬ 
gation  models  constructed  and  computed  through  the  code  correspond  to  sing¬ 
ular  integro-differentlal  equation  as  well  as  partial  differential  equation 
approximations  to  the  one-way  wave  equation.  Indeed,  this  numerical  algo¬ 
rithm  represents  one  of  the  very  few  attempts  to  compute  directly  with 
pseudo-differential  and  Fourier  Integral  operators.  For  the  two-dimensional 
case,  the  range-incrementing  procedure  Is  just  a  sequence  of  matrix  multi¬ 
plications,  and,  thus.  Ideally  suited  for  computers  which  provide  either  a 
vector  or  a  parallel  pipe  type  of  operation.  Phase  space  filtering  reduces 
both  the  size  of  the  matrix  multiplication  and  the  number  of  matrix  elements 
Initially  computed,  in  particular,  reducing  the  total  range- Incrementing 
computational  time  by  almost  an  order  of  magnitude  for  typical  model  calcu¬ 
lations  [4], 

Numerical  results  of  transmission  loss  (dB  re  1  m)  as  a  function  of 
range  (km)  for  a  number  of  model  ocean/bottom  propagation  experiments 
demonstrate  the  computational  viability  of  the  factor1zation-/path  Integra¬ 
tion-based  phase  space  marching  algorithm  [4,12].  Several  propagation 
experiments  are  summarized  In  Figures  3,  5,  and  7,  with  the  corresponding 
transmission  loss  curves  compared  with  a  reference  Fast  Field  Program  (FFP) 
algorithm  [4,12]  In  Figures  4,  6,  and  8. 


Transmission  Loss  (dB) 


Transmission 
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ATTENUATING  BOTTOM 


FT  BOTTOM 


Sound  speed  (m/»> 


Fig.  5.  Model  environment  2  and  propagation  experiment 
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7.  Model  environment  3  and  propagation  experiment. 


Fig.  8.  Transmission  loss  (dB  re  1  m)  versus  range  (km)  for  model 

environment  3  at  25  Hz.  ( - )  Real  Weyl  High  Frequency  ( - ) 
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For  400  Hz  propagation  In  the  exaggerated  double-well  Model  of  Figure 
3,  a  wide-angle  capability  well  beyond  the  ordinary  parabolic  approxlMatlon 
Is  required.  Figure  4  Illustrates  the  excellent  agreeMent  between  the  high- 
frequency  and  FFP  algorlthMS  over  ranges  on  the  order  of  500  wavelengths. 
Figure  6  Illustrates  the  cumulative  growth  of  a  phase  shift  error  at  long 
range  which  characterizes  the  breakdown  of  the  high-frequency  algorlthM  In 
the  250  Hz  propagation  In  the  rapidly  changing  shallow-water  model  of  Figure 
5.  Combining  Fourier  component,  or  wave  number,  filtering  with  the  high- 
frequency  algorithm  leads,  not  only  to  a  more  efficient  and,  thus,  faster 
algorithm,  but  also,  to  a  more  widely  applicable  numerical  scheme.  The 
filtering.  In  addition  to  removing  Fourier  components  which.  In  principle, 
make  no  significant  contribution  to  the  computed  wave  field,  eliminates 
those  unnecessary  regions  of  phase  space  where  the  small  error  In  the  high- 
frequency  symbol  approximation  can  lead.  In  a  cumulative  manner,  to  serious 
discrepancies  at  sufficiently  long  ranges.  This  Is  particularly  well 
Illustrated  In  the  60  degree  filtered  calculation  on  model  environment  2  at 
250  Hz  which  results  In  the  complete  elimination  of  the  cumulative  phase 
shift  error  (Figure  6),  greatly  extending  the  effective  computational  range. 
Sufficiently  decreasing  the  propagation  frequency  and  Increasing  the  jump 
discontinuity  in  the  sound  speed,  as  Illustrated  In  the  25  Hz  propagation  In 
the  shallow-water  model  of  Figure  7,  demonstrate  the  violation  of  energy 
conservation  Inherent  In  the  high-frequency  wave  theory  and  the  now-rapid 
decay  with  Increasing  range  of  the  corresponding  numerical  algorithm.  This 
growth  In  the  wave  field,  Illustrated  In  Figure  8,  Is  eliminated  by  the  real 
Weyl  high-frequency  algorithm  [4],  which  effectively  restores  energy 
conservation,  as  Is  also  Illustrated  In  Figure  8.  A  more  detailed 
discussion  of  these  and  other  points  Is  presented  elsewhere  [1,4,12]. 

The  speed  and  modest  storage  requirements  of  the  filtered  one-way 
algorithm  Indicate  that  range-dependent  calculations  over  extended 
environments  should  be  feasible  with  current  supercomputer  technology.  Both 
range-updating  and  the  numerical  calculation  of  the  reflected  and 
transmitted  fields  from  an  interface  should  be  possible  over  distances  on 
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the  order  of  10  wavelengths.  Preliminary  computations  with  range-dependent 
Munk-proflle  deep  ocean  environments,  Including  propagation  through  extended 
shadow  regions,  compare  well  with  adiabatic  normal -mode  calculations. 

Both  the  range-dependent  and  range- independent  Feynman/Garrod  path 
Integral  representations  can  be  computed  by  standard  Monte  Carlo  (statisti¬ 
cal  sampling)  methods  for  the  numerical  evaluation  of  multiple  Integrals 
[4].  While  numerically  calculating  Helmholtz  wave  fields  as  high  (In 
principle,  1nf1n1te)-d1mens1onal  Integrals  Is  quite  distinct  from  the  more 
traditional  finite-difference  and  finite-element  approaches,  the  Monte  Carlo 
evaluation  of  functional  Integrals  has  been  successfully  applied  In  quantum 
mechanical,  statistical  mechanical,  and  quantum  field  theoretical 
calculations  [4].  For  the  phase  space  representations  of  Eqs.  (13)  and  (16) 
in  two  dimensions  (n  *  2),  the  modeling  of  realistic  propagation  experiments 
can  Involve  the  computation  of  thousand-dimensional  oscillatory  Integrals. 
Correlated-sampllng  variance  reduction  techniques  can  dramatically  Improve 
the  speed  and  accuracy  of  the  algorithm  [4].  Generally  speaking,  a  large 
parallel  processing  capability  should  have  a  very  favorable  Impact  on  the 
numerical  computation  of  path  Integrals  [4]. 

IV.  INVERSE  FORMULATION.  The  phase  space-based  construction  of  the 


square  root  Helmholtz  operator  provides  the  basis  for  a  formulation  of  the 
Inverse  algorithms  mentioned  In  the  Introduction.  Mathematically,  the 
refractive  Index  field  (or  Its  square)  Is  reconstructed  from  the  full-space 
Ke1mho1t2  Green's  function  G  through  the  relationship 

B(xt,xM  *  (21/C )  11m  to*G(x,xJO,xi)).  (26) 

-t  -i  x — >0  x  ~z 

The  symbol  Og(£,<})  Is  then  constructed  through  an  Inverse  Fourier  transform 
of  the  kernel  function  B(xt,x^)  and  subsequently  yields  the  refractive  Index 

field  upon  a  direct  application  of  the  Weyl  composition  equation  (3)  forlpl* 

0.  In  the  homogeneous  medium  limit,  the  direct  evaluation  of  the  composite 

o 

symbol  reduces  to  the  square  of  the  symbol,  Qg2(£,<i)  *  ^(£,5). 

The  Inverse  algorithm  proceeds  around  the  correspondence  diagram  (pictorial 
summary)  In  a  counterclockwise  fashion.  The  direct  propagation  algorithm 
requires  the  Inversion  of  Eq. (3 )  while  the  Inverse  propagation  algorithm 
only  requires  a  direct  computation  of  Eq.(3).  Thus  the  direct  propagation 
problem  has  been  transformed  Into  an  "Inverse"  problem  while  the  wave  field 
Inversion  problem  has  been  reformulated,  In  an  appropriate  sense,  as  a 
direct  calculation. 

The  factorization  algorithm  exactly  Inverts  the  Inherently  nonlinear 
relationship  between  the  wave  field  data  and  the  refractive  Index  field  as 
reflected  In  the  Llppmann-Schwlnger  equation  for  the  propagator  [3].  Most 
importantly.  It  Is  a  multidimensional  formulation.  For  the  "physical 
experiment,"  a  point  source  Is  Introduced  Into  the  medium  defining  the 
Initial-value  (x  ■  0)  plane.  The  second  derivative  with  respect  to  the 
range  of  the  wave  field  Is  then  determined  as  a  function  of  the  point  source 
and  receiver  positions.  Collecting  the  data  on  the  Initial -value  plane 
would  most  probably  limit  the  application  of  the  algorithm  to  specific  types 
of  bore-hole  experiments.  Moreover,  mathematically,  the  Inversion  requires 
the  evaluation  of  singular  Integrals  (generalized  functions).  Collecting 
data  on  a  downfleld  plane  (x  >  0)  leads  to  a  transmission  experiment  similar 
to  the  oceanic  sound  speed  profile  Inversion  method  of  DeSanto  [3].  The 
downfleld  wave  field  provides  for  an  appropriate  analytic  continuation  In 
the  factorization  algorithm  and  connects  the  analysis  with  the  Inverse 
diffraction  problem  [3]. 

The  transmission,  or  propagation,  formulation  Is  analogous  to 
tomography.  The  reference  wave  number  In  the  factorization  analysis 
corresponds  to  2«/(Planck's  constant)  as  opposed  to  Its  square  playing  the 
role  of  an  energy.  The  source  generation  and  data  collection  over  parallel 
planes  then  naturally  correspond  to  the  multidirectional  Insonl fy Ing  plane 
waves  and  subsequent  angular  data  collection  of  fixed-energy  (frequency) 
diffraction  tomography  [3].  For  range- dependent  environments,  the  Inclusion 
of  backscatter  effects,  even  In  an  approximate  manner,  would  then  provide 
the  basis  for  a  generalized  acoustic  tomography,  extending  the  diffraction 
algorithms  based  on  the  Born,  Rytov,  or  distorted-wave  Born  approximations 
[3].  The  nonlinear  factorization  and  subsequent  weak -back scatter 
perturbation  theory  would  extend  the  linearized  weak -scattering  treatments 
Into  the  nonlinear  regime.  This  can  be  attempted  In  two  ways.  Formal  field 
splitting  analysis  provides  the  basis  for  a  weak -back scatter  perturbation 
theory  within  the  framework  of  Invariant  Imbedding  [2-4].  The  arbitrary- 


dimensional  nature  of  the  factorization  analysis  In  conjunction  with 
mathematical  Imbedding  concepts  provides  the  basis  for  a  spatial-dimensional 
perturbation  theory  [2-4].  This  essentially  Involves  treating  the  spatial 
dimension  of  both  the  Helmholt2  operator,  in  general,  and  the  refractive 
Index  field.  In  particular,  as  a  variable  and  subsequently  studying  the 
structure  of  the  resulting  family  of  systems  Indexed  In  this  manner.  For 
the  case  of  two  (different)  transversely  Inhomogeneous  half-spaces 
separated  by  a  planar  Interface,  an  Inverse  algorithm  can  be  Initially  based 
on  the  composition  equation  (9). 

For  a  transversely  1 nnomogeneou s  environment,  the  factorization 
Inversion  model  Invites  comparison  with  "effective  one-dimensional" 
stratified  environmental  models  such  as  that  of  Stickler  and  Delft  [4].  In 
both  models,  the  location  of  the  field  source  (finite)  and  the  data  measure¬ 
ments  Is  within  the  scattering  region.  Most  Importantly,  the  factorization 
method  Is  a  direct  Inversion  of  an  arbitrary-dimensional  propagation 
equation  which  requires  less  symmetry  than  those  models  (l.e., 
Stickler-Delft)  reducible  to  the  standard  one-dimensional  formulation  of 
Delft-Trubowltz  [4]  or  Gel fand-LevI tan  [4].  Thus  for  example.  In  a  general 
n-dlmenslonal  Cartesian  formulation,  the  refractive  Index  field  can  be  a 
function  of  as  many  as  (n-1)  coordinates  In  the  factorization  model,  while  a 
function  of  only  one  coordinate  In  an  "effective  one- dimensional "  model. 

The  experiment  envisioned  and  the  distinguished  direction  differ  In  the  two 
models.  In  the  transversely  Inhomogeneous  environment,  the  direction  In 
which  there  Is  medium  homogeneity  Is  distinguished,  while  In  the  "effective 
one-dimensional"  model,  the  one  direction  In  which  there  Is  medium  Inhomo¬ 
geneity  Is,  In  effect,  distinguished.  Data,  In  both  cases.  Is  collected 
perpendicular  to  the  distinguished  direction.  The  Stickler-Delft  model  Is 
essentially  a  one-dimensional  scattering  experiment  with  the  surface  data. 

In  effect,  reflection  coefficient  data.  Thus  unlike  the  transmission 
experiment,  which  extensively  samples  the  region  of  Inhomogeneity,  In  the 
factorization  model,  the  Stickler-Delft  analysis  does  not  account  for  the 
presence  of  "trapped  modes"  [4].  The  formal  Inclusion  of  a  specific 
pressure-release  surface  within  the  pseudo-differential  operator  framework 
would  allow  for  a  stratified  environmental  model  and  the  subsequent 
quantitative  comparison  with  the  Stickler-Delft  model. 

For  applied  Inverse  problems,  approximate  Inversions  may  prove 
adequate.  Approximate  Inversion  algorithms  follow  readily  from  the  pertur- 

batlve  treatments  of  the  Ueyl  composition  equation.  K  (3)  Is  related  to 
QgfO,^)  In  a  quadratic  fashion  and  through  a  linear  Integral  relationship, 

respectively.  In  the  high-frequency  (£ — >•)  and  weak-lnhomogenelty  (Born) 
limits.  In  particular,  the  high-frequency  algorithm  Is  based  upon  choosing. 
In  practice,  a  I p  I  such  that  the  symbol  approaches  Its  asymptotic  form, 

2  2  1/2 

(K  (<|)  -  £  )  •  The  approach  to  the  asymptotic  regime  In  phase 

space  Is  governed  both  by  the  magnitude  of  KZ(<|)  -  p2  (large)  and  the 
variation  of  the  refractive  Index  field  on  the  wavelength  scale  (small). 
Figure  9  Illustrates  the  high-frequency  inversion  for  the  case  of  a 
quadratic  medium.  Applying  the  full  composition  equation  for  the  Inversion 
would  result  In  a  linear  function  In  X  for  the  real  part  and  an  Imaginary 
part  which  Is  Identically  zero.  Finally,  weighted  Hilbert  space  methods  for 
Incorporating  prior  estimates  appear  to  be  applicable  to  the  Fourier-based 
factorization  approach  [4]. 
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